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ABSTRACT 

A revision of Stodoikiewicz’s Monte Carlo code is used to simulate the evolution of million 
body star clusters. The new method treats each superstar as a single star and follows the 
evolution and motion of all individual stellar objects. A survey of the evolution of A^-body 
systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. 
The process of energy generation is realized by means of appropriately modified versions of 
Spitzer’s and Mikkola’s formulae for the interaction cross section between binaries and field 
stars and binaries themselves. The results presented are in good agreement with theoretical 
expectations and the results of other methods. During the evolution, the initial mass function 
(IMF) changes significantly. The local mass function (LMF) around the half-mass radius 
closely resembles the actual global mass function (GMF). At the late stages of evolution the 
mass of the evolved stars inside the core can be as high as 97% of the total mass in this 
region. For the whole system, the evolved stars can compose up to 67% of the total mass. The 
evolution of cluster anisotropy strongly depends on initial cluster concentration, IMF and the 
strength of the tidal field. The results presented are the first step in the direction of simulating 
the evolution of real globular clusters by means of the Monte Carlo method. 

Key words: gravitation - methods: numerical - celestial mechanics, stellar dynamics - glob¬ 
ular clusters: general 


1 INTRODUCTION. 

Dynamical modeling of real, large stelar systems, like globular 
clusters or galactic nuclei, and understanding their evolution still 
is a great challenge both for the theory and hardware/software. Ba¬ 
sically, there are two approaches. Direct A-body method, which 
requires extremely large hardware requirements and very sophis¬ 
ticated software, and statistical modeling based on the Fokker- 
Planck and other approximations, which suffers from the poor un¬ 
derstanding of the validity of assumptions. 

On the side of A-body simulations recent years brought a 
significant progress in both hardware and software. First, parallel 
(even massively parallel) computing has opened a route to gain 
performance at relatively low cost and little technological advance¬ 
ment. Secondly, the already successful GRAPE special purpose 
computers (Sugimoto et al. 1990, Makino et al. 1997, Makino & 
Taiji 1998, Makino 2002, Makino & Hut 2003) have been de¬ 
veloped in their present generation (GRAPE-6) and aim at the 
100 Tflops speed. NBODY6-b+ (Spurzem 1999), the successor of 
Aarseth’s NBODY6 code has been ported on massively parallel, 
general purpose computers (CRAY T3E and PC clusters) and was 
used for astrophysical problems related to the dissolution of glob- 
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ular clusters (Baumgardt 2001, Baumgardt, Hut & Heggie. 2002) 
and the decay of massive black hole binaries in galactic nuclei after 
a merger (Milosavljevic & Merritt 2001, Hemsendorf, Sigurdsson 
& Spurzem 2002). Despite such progresses in hardware and soft¬ 
ware there is still impossible to model directly evolution of real 
globular cluster (A ~ 10®) and galactic nuclei (A ~ 10®). Recent 
work by Baumgardt et al (2002) and Baumgardt & Makino (2002) 
has pushed the limits of present direct modeling to about 10® using 
both NBODY6-H- and the GRAPE-6 special purpose hardware. 

On the side of Fokker-Planck method with finite differences 
and an anisotropic gaseous model, recent years brought great im¬ 
provements. Models can be used now to simulate more realistic 
stellar systems. They can tackle: anisotropy, rotation, a tidal bound¬ 
ary, tidal shocking by galactic disk and bulge, mass spectrum, 
stellar evolution and dynamical and primordial binaries (Louis & 
Spurzem 1991, Giersz & Spurzem 1994, Takahashi 1995, 1996, 
1997, Spurzem 1996, Takahashi & Portegies Zwart 1998, 2000 
hereafter TPZ, Drukier et al 1999, Einsel & Spurzem 1999, Taka¬ 
hashi & Lee 2000, Kim et al. 2002, Kim, Lee & Spurzem 2004, 
Fiestas, Spurzem & Kim 2005). Unfortunately, the Fokker-Planck 
approach suffers, among other things, from the uncertainty of dif¬ 
ferential cross-sections of many processes which are important 
during cluster evolution. It can not supply detailed information 
about the formation and movement of all objects present in clus- 
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ters. Additionally, the anisotropic gaseous model assumes a certain 
form of heat conductivity and closure relations between the third 
order moments. 

There is an elegant alternative to generate models of star clus¬ 
ters, which can correctly reproduce the stochastic features of real 
star clusters, but without really integrating all orbits directly as in 
an Al-body simulation. They rely on the Fokker-Planck approxi¬ 
mation and (hitherto) spherical symmetry, but their data structure 
is very similar to an A^-body model. These so-called Monte Carlo 
models were first presented by Henon (1971, 1975), Spitzer (1975) 
and later improved by Stodolkiewicz (1982, 1985, 1986) and in fur¬ 
ther work by Giersz (1997, 1998, 2001), and recently by Rasio and 
his collaborators (Joshi, Rasio & Portegies Zwart 2000, Watters, 
Joshi & Rasio 2000, Joshi, Nave & Rasio 2001, Fregeau et al. 2003, 
Giirkan, Freitag & Rasio 2004), and Freitag (Freitag 2000, Freitag 
& Benz 2001, 2002). The basic idea is to have pseudo-particles, 
which orbital parameters are given in a smooth, self-consistent po¬ 
tential. Flowever, their orbital motion is not explicitly followed; to 
model interactions with other particles like two-body relaxation by 
distant encounters or strong interactions between binaries and other 
objects, a position of the particle in its orbit, and further free param¬ 
eters of the individual encounter, are picked from an appropriate 
distribution by using random numbers. The Monte Carlo scheme 
takes full advantage of the established physical knowledge about 
the secular evolution of (spherical) star clusters as inferred from 
continuum model simulations. Additionally, it describes in a proper 
way the graininess of the gravitational field and the stochasticity of 
real A-body systems and provides, in a manner as detailed as in di¬ 
rect A-body simulations, information about the movement of any 
objects in the system. This does not include any additional physi¬ 
cal approximations or assumptions which are common in Fokker- 
Planck and gas models (e.g. for conductivity). Because of this, the 
Monte Carlo scheme can be regarded as a method which lies be¬ 
tween direct A-body and Fokker-Planck models and combines 
most of their advantages. A hybrid variant of the Monte Carlo 
technique combined with a gaseous model has been proposed by 
Spurzem & Giersz (1996), and applied to systems with a large num¬ 
ber of primordial binaries by Giersz & Spurzem (2000) and Giersz 
& Spurzem (2003) and for tidally limited systems by Spurzem et 
al (2005). The hybrid method uses a Monte Carlo model for bi¬ 
naries or any other object for which a statistical description, as 
used by the gaseous model, is not appropriate due to small num¬ 
bers of objects or unknown analytic cross sections for interaction 
processes. The method is particularly useful for investigating the 
evolution of large stellar systems with a realistic fraction of pri¬ 
mordial binaries, but could also be used in the future to include, for 
example, the huild up of massive stars and blue stragglers by stellar 
collisions. 

Very detailed observations of globular clusters which are 
available at present or will become available in future (e.g. the Hub¬ 
ble Space Telescope and the new very large terrestrial telescopes) 
have extended and will extend our knowledge about their stellar 
content, internal dynamics and the influence of the environment on 
cluster evolution (Janes 1991, Djorgovski & Meylan 1993, Smith 
& Brodie 1993, Hut & Makino 1996, Meylan & Heggie 1997). 
This data covers luminosity functions and derived mass functions, 
color-magnitude diagrams, population and kinematical analysis, 
including binaries and compact stellar evolution remnants, detailed 
two-dimensional proper motion and radial velocity data, and tidal 
tails spanning over arcs several degrees wide (Koch et al. 2004). 
A wealth of information on ’’peculiar” objects in globular clusters 
(blue stragglers. X-ray sources (high- and low-luminosity), mil¬ 


lisecond pulsars, etc.) suggests a very close interplay between stel¬ 
lar evolution, binary evolution and dynamical interactions. This in¬ 
terplay is far from being understood. Moreover, recent observations 
suggest that the primordial binary fraction in a globular cluster can 
be as high as 15% - 38% (Rubenstein & Bailyn 1997), and possible 
existence of intermediate-mass black holes in the centers of some 
globular clusters (Gebhardt et al 2000). With all these observational 
data such as King at al. (1998), Piotto & Zoccali (1999), Ruben¬ 
stein & Bailyn (1999), Ibata et al. (1999), Piotto et al. (1999), Grill- 
mair et al. (1999), Shara et al. (1998), Odenkirchen et al. (2001), 
Hansen et al. (2002), Richer et al. (2002) (to mention only a few 
papers), easily reproducible reliable modelling becomes more im¬ 
portant than before. Monte Carlo codes provide all the necessary 
flexibility to disentangle the mutual interactions between all physi¬ 
cal processes which are important during globular cluster evolution 
and to perform in a reasonable time, detailed simulations of real 
globular clusters. 

The ultimate aim of the project described here is to build a 
Monte Carlo code (in the framework of the MODEST international 
collaborations, http:/www.manybody.org/inodest), which will be 
able to simulate the evolution of real globular clusters, as closely as 
possible. In this paper (the third in the series) the Monte Carlo code 
(which was discussed in detail in Papers I and II) is used to simulate 
the evolution of stellar systems, which consist of comparable num¬ 
ber of stars and have mass comparable to real globular clusters. The 
results of simulations will be compared with those of Chemoff & 
Weinberg (1990, hereafter CW), Vesperini c& Heggie (1997, here¬ 
after VH), Aarseth & Heggie (1998, hereafter AH), Baumgardt & 
Makino (2002) and Tamers et al (2005). 

The plan of the paper is as follows. In Section 2, a short de¬ 
scription of processes implemented into the Monte Carlo code will 
be presented. In Section 3, the initial conditions will be discussed 
and results of the simulation will be shown. And finally, in Section 
4 the conclusions will be presented. 


2 MONTE CARLO METHOD. 

The Monte Carlo method can be regarded as a statistical way of 
solving the Fokker-Planck equation. Its implementation presented 
in Giersz (1998, hereafter Paper I) and Giersz (2001, hereafter Pa¬ 
per II) is based on the orbit-averaged Monte Carlo method devel¬ 
oped in the early seventies by Henon (1971) and then substantially 
improved by Stodolkiewicz (1986, and references therein), and re¬ 
cently by Giersz (Paper I and II), by Freitag (Freitag & Benz 2001, 
2002) and by Rasio’s group (Joshi et al 2000, 2001, Fregeau et 
al 2003, 2004, Fregeau, Chatterjee & Rasio 2005, Freitag, Giirkan 
& Rasio 2005, Freitag, Rasio & Baumgardt 2005). The code is de¬ 
scribed in detail in Paper I, which deals with simulations of iso¬ 
lated, single-mass systems and in Paper II in which additional 
physical processes were included: multi-mass systems and stellar 
evolution, mass loss through the tidal boundary, formation of three- 
body binaries and their subsequent interactions with field stars, and 
interactions between binaries. 

For the sake of completeness, a short description of the main 
ingredients of the Monte Carlo model will follow. For more detailed 
discussion of the code readers should refer to Paper I and II. 

2.1 Mass Spectrum and Stellar Evolution 

To facilitate comparison with results of other numerical simulations 
of globular cluster evolution (e.g. CW, Fukushige & Heggie 1995, 
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Giersz & Heggie 1997, VH, AH, TPZ, JNR) and for simplicity a 
single power-low approximation of IMF was used. 

N{m)dm = Cm~°‘dm, rrimin ^ m ^ rtimax, (1) 

where C and a are constants. The rrimin, minimum stellar mass 
and rtimax, maximum stellar mass were chosen to O.IM© and 
I 5 .OM 0 , respectively. However, it should be noted that observa¬ 
tions give more and more evidence that the mass function in glob¬ 
ular/open clusters and for field stars as well is not a simple power- 
law, but it could be rather approximated by a composite power-law 
(e.g. Kroupa, Tout & Gilmore 1993, Kroupa 2002). To describe 
the mass loss due to stellar evolution the same simplified stellar 
evolution model as adopted by CW was used. More sophisticated 
models of stellar and binary evolution based on fitting analytical 
formulae to the evolution tracks for stars and binaries were devel¬ 
oped by Portegies Zwart & Verbunt (1996), Hurley, Pols & Tout 
(2000), Hurley, Tout & Pols (2003) and Belczynski, Kalogera & 
Bulik (2002). These models, among others, were used in population 
synthesis codes (e.g. Belczynski & Taam 2004), and simulations of 
the evolution of real star cluster M67 (e.g. Hurley et al 2005). 

The initial masses of stars are generated from the continuous 
distribution given in equation (1). This ensure a natural spread in 
their lifetimes. To ensures that the cluster remains close to virial 
equilibrium during rapid mass loss due to stellar evolution, special 
care is taken, that no more than 3% of the total cluster mass is lost 
during one overall time-step. 


2.2 Tidal Stripping 

For tidally limited star cluster the rate of mass loss is much higher 
than for isolated systems. This is connected with the tidal stripping 
of mass across the system boundary by the gravitational field of a 
parent galaxy. For isolated system mass loss is attributed mainly to 
rare strong interactions in the dense, inner part of the system. In the 
present Monte Carlo code a mixed criterion is used to identify es¬ 
capers: a combination of apocentre and energy-based criteria (see 
Paper II). This means that no potential escapers are kept in the sys¬ 
tem. Stars regarded as escapers are instantaneously lost from the 
system. This is in contrast to A^-body simulations where stars need 
time proportional to the dynamical time to be removed from the 
system. Recently, Baumgardt (2000) showed that stars with energy 
greater than the tidal cutoff energy Et = —GM/rt, in the course 
of escape, can again become bound to the system, because of dis¬ 
tant interactions with field stars. A similar conclusion was reached 
by Fukushige & Heggie (2000), who showed that the escape time- 
scale can be very long (much larger than the half-mass relaxation 
time). Additionally, during the final stages of cluster evolution or 
for clusters with initially low central concentration, the mass loss 
across the tidal boundary can become unstable when too many stars 
are removed from the system at the same time. To properly follow 
these stages of evolution the time-step has to be decreased to force 
smaller mass loss. 


2.3 Three-Body Binaries and their interactions 

In the present Monte Carlo code (as it was described in Paper I) all 
stellar objects, including binaries, are treated as single superstars. 
This allows to introduce into the code, in a simple and accurate 
way, processes of stochastic formation of binaries and their sub¬ 
sequent stochastic interactions with field stars and other binaries. 
The procedure of binary formation in three-body interactions is in 


great detail described in Paper II. It relies on the observation that 
the probability that the masses of the three stars involved in the 
interaction are mi, m 2 and m 3 is nin 2 n 3 /n® (where ni, 712 , ns 
and n are number densities of three interacting stars and the total 
number density, respectively) and the rate of binary formation is 
proportional to n\n 2 ns. These considerations lead to the formula 
for the probability of binary formation which depends only on the 
local total number density instead of local number densities of each 
mass involved in the interaction (see Equation 7 in Paper II). This 
procedure substantially reduces fluctuations in the binary formation 
rate. The determination of the local number density in the Monte 
Carlo code, particularly for multi-mass systems, is a very delicate 
and difficult matter (Paper I). 

The energy generation in binary - field star interactions is 
computed according to the (appropriately modified for multi-mass 
case) semi-empirical formula of Spitzer (Spitzer 1987). Unfor¬ 
tunately, for multi-component systems there is no simple semi- 
empirical formula which can fit all numerical data (Heggie, Hut & 
McMillan 1996). The total probability of the interaction between 
a binary of binding energy Et, consisting of mass mi and m 2 and 
a field star of mass m 3 was deduced from Equations (6-23), (fi¬ 
ll) and (6-12) presented in Spitzer (1987) and results of Heggie 
(1975) with an appropriately adjusted coefficient (the total proba¬ 
bility computed is such a way gives correct value for the equal- 
mass case). 


^sLi-KAsG'^mlmlCrnUsn 

^3bf = - — - -At, (2) 

8v 3^mi2 yTniyTnlcrAi, 

where mi 2 = mi-|-m 2 , mi 23 =mi 2 -|-m 3 , m^ is the mean stellar 
mass in a zone, and a is the one dimensional velocity dispersion. 
This procedure is, of course, oversimplified and in some situations 
can not give correct results, for example, when a field star is very 
light compared to the mass of the binary components. 

The implementation of interactions between binaries in the 
code is based on the method described by Stodolkiewicz (1986) 
and then improved in Paper II. Only strong interactions are consid¬ 
ered, and only two types of outcomes are permitted: heavier binary 
and two single stars, or two binaries in a hyperbolic relative or¬ 
bit. Stable three-body configurations are not allowed. In the case 
when binaries are only formed in dynamical processes, and only a 
few binaries are present, at any time, in the system, it is very diffi¬ 
cult to use the binary density (Giersz & Spurzem 2000) to calculate 
the probability of a binary-binary interaction. Given binary can hit 
only binaries (regarded as targets) whose actual distance from the 
cluster center lies between its pericenter and apocenter. This leads 
to the following formula for the binary-binary interaction probabil¬ 
ity. 


wp 

^ Iml 2r2r^^’ 


(3) 


where w is the relative velocity between interacting binaries, Iw^l, 
r, T and p are the radial velocity, radial position, orbital period of 
a binary in the system and impact parameter (equal to 2.5 times the 
semi-major axis of the softer binary), respectively. The outcome 
of the interaction is as follows (see details in Stodolkiewicz 1986, 
Paper II and Mikkola 1983, 1984): 


• Two binaries in a hyperbolic relative orbit (12% of all interac¬ 
tions). Energy generation is equal to 0.4 times the binding energy 
of the softer binary. 

• One binary and two escapers ( 88 % of all interactions). The 
softer binary is destroyed and the harder binary increases its bind- 
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Table 1. Models “ 


Model 

Wo 

a 

Adj’ 

Rg 


^scale 

F 

W3235 

3 

2.35 

319305 

2.00 

3.1311 

79797 

4.62 

W335 

3 

3.5 

166577 

2.77 

3.1311 

61418 

3.34 

W5235'’ 

5 

2.35 

319305 

2.00 

4.3576 

48602 

4.62 

W7235® 

7 

2.35 

319305 

2.00 

6.9752 

23999 

4.62 


“ Mj- is the total cluster mass in solar mass, rt/^ is the tidal 
radius in Monte Carlo units, tg^ale is the time scaling factor to 
scale simulation time to physical time in 10® yrs, Rq is the 
distance from the Galactic Center in kpc, F in 10^ (Chemoff 
& Weinberg 1990). The first entry after W describes the King model 
and the following numbers the mass function power-law index. 

^ two models with different initial random seed were simulated. 


ing energy by the amount equal to 0.516 times the sum of binding 
energies of binaries. 

The procedures described above are very uncertain in regard to 
the amount of energy generation. To cure this problem, it is planned 
in the near future to introduce into the code numerical procedures 
(based on Aarseth’s NBODY 6 code), which can numerically in¬ 
tegrate the motion of three- and four-body subsystems. A simi¬ 
lar procedure was introduced, with success, into the Hybrid code 
(Giersz & Spurzem 2000, 2003). 


3 RESULTS 

In Paper II a survey of the evolution of A^-body systems influenced 
by the tidal held of a parent galaxy, by stellar evolution and by en¬ 
ergy generation in interactions between binaries and field stars and 
binaries themselves was discussed. Here, the first Monte Carlo sim¬ 
ulations of million body systems are presented. The results will be 
discussed from the point of view of system parameters which can 
be verified by observational data. Among others, evolution of LMF 
and GMF, anisotropy, total mass, mean mass, mass segregation and 
density profiles will be presented. 


3.1 Initial Models 

The initial conditions were chosen in a similar way as in the “col¬ 
laborative experiment” (Heggie et al 1999) and Paper II. The posi¬ 
tions and velocities of all stars were drawn from the King model 
(King 1966). All models have the same total number of stars 
Nt = 1, 000, 000 and the same tidal radius rt = 33.57 pc. Masses 
of stars are drawn from the power-law mass function according to 
Equation ( 1 ). The minimum mass was chosen to be O.IMq and the 
maximum mass 15M0. Two different values of the power-law in¬ 
dex were chosen: a = 2.35 and 3.5. The sets of initial King models 
were characterized by Wo = 3, 5 and 7. All models are listed in 
Table 1. 

F is the parameter (introduced by CW), which is proportional 
to the half-mass relaxation time. 

^ _ Mt Rg 220 kms~^ 1 
Mq kpc Vg InN ’ 

where Rq is the distance of the globular cluster to the Galactic 
Center, and Vg is the circular speed of the cluster around the Galaxy. 


In order to scale the Monte Carlo time to physical units the follow¬ 
ing formula is used 


iscale 

Kfiyrs 


14.91 


\ 1-5 

n \ Nt 

rtM ) V MtIuGNt) ’ 


(5) 


where rt^ is the tidal radius in Monte Carlo units, M is the total 
cluster mass and rt is the tidal radius, both are in solar mass and 
pc, respectively and 7 is equal to 0.11 as for single-mass systems 
(Giersz & Heggie 1994). However, for multi-mass systems it could 
be much smaller (Giersz & Heggie 1996, 1997). The initial model 
is not in virial equilibrium, because of statistical noise and because 
of masses, which are assigned independently from the positions 
and velocities. Therefore, the model has to be initially rescaled to 
fulfill the virial equilibrium condition. During all the simulations 
the virial ratio is kept within < 2 % of the equilibrium value (for 
the worst case within < 5%). Standard A-body units (Heggie & 
Mathieu 1986), in which the total mass M — 1,G = 1 and the ini¬ 
tial total energy of the cluster is equal to —1/4, have been adopted 
for all runs. Monte Carlo time is equal to A-body time divided by 
Nt/ItiGNt)- In the course of evolution, when the cluster looses 
mass, the tidal radius changes according to rt oc 

Finally, a few words about the efficiency of the Monte Carlo 
code presented here. The simulation of a million body systems took 
about two/three weeks (depending on an initial model) on a Pen¬ 
tium 4, 2 GHz PC. This is still much faster than the biggest direct 
A-body simulations performed on GRAPE-4 (Teraflop special- 
purpose hardware). Nevertheless, the speed of the code is not high 
enough to perform large scale survey simulations in a reasonable 
time. It is clear, that to simulate the evolution of real globular clus¬ 
ters with all physical processes in operation, a substantial speed-up 
of the code is needed. This can be done, either by parallelizing the 
code (in a similar way to INR), introducing a more efficient way 
of determining the new positions of superstars, introducing super- 
stars which contain a different number of stars, or by using a hybrid 
code (as was done by Giersz & Spurzem 2000, 2003). However, the 
Monte Carlo code has presently a great relative advantage over the 
A-body code for simulations with a large number of primordial 
binaries. Primordial binaries substantially downgrade the perfor¬ 
mance of A-body codes on supercomputers or on special-purpose 
hardware. They can be introduced into the Monte Carlo code in a 
natural way, practically without a substantial loss of performance 
(Giersz & Spurzem 2000, 2003). 


3.2 Global evolution. 

The rate of cluster evolution and indirectly the strength of the tidal 
field is described by the parameter F (Equation. 4). Generally, the 
greater the value of F, the longer the relaxation time and the slower 
the cluster evolution. The same relation is true also for the distance 
of a globular cluster to the Galactic Center, Rg (keeping the con¬ 
stant total cluster mass Mt and number of stars At). The greater 
the F, the larger the Rg- As one can see from the Table 1 dis¬ 
tance Rg is rather small This means that the galactic tidal field is 
very strong and globular cluster evolution is very fast. These ini¬ 
tial conditions were chosen mainly to speed up the simulations and 
discussed models similar to models considered by VH. 

Figure 1 presents the evolution of the total mass for all dis¬ 
cussed models. Three phases of evolution are clearly visible. The 
first phase is connected with the violent mass loss due to stellar 
evolution of the most massive stars. Then, there is the long phase 
connected with gradual mass loss due to tidal stripping. And fi- 
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Figure 1. Evolution of the total mass for models W7235, W5235, W3235 Figure 2. Mass loss due to stellar evolution (evo) and tidal stripping (esc) 

and W335. for models W7235, W5235, W3235 and W335 


nally, there is the phase connected with the tidal disruption of a 
cluster. The last phase is well visible only in model W7235, which 
evolution is followed long enough. In model W3235 on can see 
a phase when the total cluster mass is nearly constant and very 
small (about a few thousands of the initial mass). This is an arti¬ 
ficial effect connected with the fact that the Monte Carlo code can 
not properly deal with systems with the total binding energy close 
to zero and consist of only a few dozen particles. It is interesting to 
note, that the second phase (gradual tidal stripping) can be divided 
into two subphases, clearly separated by the core collapse time. 
For the post collapse phase the rate of mass loss is slightly steeper 
than in the collapse phase. This can be connected with the fact that 
hard interactions, which involve binaries, start to operate introduc¬ 
ing an additional mechanism of the mass loss (escaping stars and 
binaries). Models with the same mass-function (a ^ —2.35) and 
with different Wo show very similar evolution during the phase of 
tidal stripping (except W3235). It seems that the initial mass loss 
due to stellar evolution is not sufficiently strong to substantially 
change the initial cluster structure. Model W3235 is not initially 
concentrated strongly enough to survive without substantial struc¬ 
ture changes, the violent mass loss event. Figure 1 presents qual¬ 
itatively very similar features as figures shown by VH, Paper II, 
Baumgardt & Makino (2002) and Tamers et al (2005), Figure I, 
Figures 1-3, Figure I and Figure 2, respectively. 

The total mass loss consists of mass loss due to stellar evolu¬ 
tion (evo) and mass loss connected with tidal stripping (esc). Re¬ 
sults presented in Figure 2 are in a good qualitative agreement with 
the A^-body results of VH (Figure 2), taking into account that a 
adopted by VH is —2.5 instead —2.35. According to Equation 11 
of VH, the evolution of the total mass of the cluster can be fitted 
by a simple analytical expression, which consists of two terms: evo 
and esc. The amount of mass loss because of esc can be fitted (as it 
was shown in VH) by a straight line with a coefficient (3. The val¬ 
ues of 13 listed by VH are 0.828 and 0.790 for models Wo = 7, and 
Wo = 5, respectively. The equivalent values for the Monte Carlo 
models are 0.999 and 0.943. The bigger values of (3 for the Monte 
Carlo models can be attributed mainly to the different: escape cri¬ 
terion, power-low index of the initial mass function and slight in¬ 
crease of the rate of mass loss due to binary activities (see Figure 1). 
As it was stated in Section 2.2, the escape criterion adopted in this 
paper (together with the instantaneous removal of escapers) leads to 
a higher escape rate than in A^-body simulations. Additionally, less 


steep initial mass function causes stronger mass loss due to stellar 
evolution and consequently leads to more escapers than the steeper 
mass function adopted by VH. Mass loss connected with stellar 
evolution dominates the initial phase of cluster evolution. Then the 
rate of stellar evolution substantially slows down and escape due 
to tidal stripping takes over. During this phase of evolution the rate 
of mass loss due to stellar evolution is nearly constant, and higher 
for shallower mass functions (as expected). Energy carried away 
by stellar evolution events dominates the energy loss due to tidal 
stripping, even though the tidal mass loss is higher. 

It is worth to note, that all models except W3235 show clear 
gravothermal oscillations in the post collapse phase. The oscilla¬ 
tions become less and less visible when the cluster is approaching 
the final stages of dissolution. 

Generally, agreement throughout the evolution between the re¬ 
sults presented here and those by AH, VH, Baumgardt c& Makino 
(2002) and Tamers et al (2005) for V-body models, TPZ for 2-D 
Eokker-Planck models and by JNR and Paper II for Monte Carlo 
models is rather good. In all cases, the qualitative behaviour is very 
similar. Nevertheless, the treatment of escapers for a tidally lim¬ 
ited system proposed by Spurzem et al (2005) can be implemented 
to properly ando more accurately follow the mass loss due to tidal 
stripping and the final stages of the cluster evolution, just before 
the dissolution. 

At the end of this section the evolution of the density profiles 
for model W7235 will be discussed (see Eigure 3). The very large 
number of particles in models allows substantial reduction in sta¬ 
tistical fluctuations and computation (with high precision) of the 
density profiles in the course of cluster evolution. From the theo¬ 
retical considerations it follows that, during the collapse phase the 
nearly isothermal core is building up. At the time of core bounce 
(about 3.5 Gyr) the density profile reaches the theoretically pre¬ 
dicted power-low index equal to —2.2. Then this index is preserved 
during the further post-collapse evolution. For models W7235, 
W5235 and W335 the evolution of density profiles is practically 
the same. Only for model W3235, which undergoes strong mass 
loss due to stellar evolution and tidal stripping the density profiles 
do not show any signs of power-low core development. The initial 
shape of the profile is practically preserved (taking into account the 
reduction of central density and tidal radius). 





6 Mirek Giersz 



log(r) 


Time (Myr) 


Figure 3. Evolution of the density profiles for model W7235. The times (in 
Gyr) at which density profiles are plotted are: 0, 2, 3, 3.5, 5, 6 and 7. The 
theoretically predicted density profile in the core plotted as a straight line 
labeled by -2.2. 


3.3 Anisotropy evolution. 

The degree of velocity anisotropy is measured by a quantity 


A = 2 



( 6 ) 


where ar and at are the radial and tangential one-dimensional ve¬ 
locity dispersions, respectively. For isotropic systems A is equal to 
zero. Systems preferentially populated by radial orbits have A pos¬ 
itive and systems preferentially populated by tangential orbits have 
A negative. The evolution of anisotropy is presented in Figures 4 to 
6 for models W3235, W5235 and W7235, respectively. As can be 
seen in Figure 4, for systems which can not survive the violent ini¬ 
tial mass loss without strong structural changes (see discussion in 
the previous section), anisotropy stays close to zero with very large 
fluctuations. The system does not live long enough and the tidal 
stripping is strong enough to prevent the developing of substantial 
positive or negative anisotropy. For model W5235 the cluster is not 
initially very strongly concentrated, so the mass loss due to stellar 
evolution and consequently by tidal stripping prevents the devel¬ 
opment of positive anisotropy in the outer and middle parts of the 
system (see Figure 5). From the very beginning anisotropy for these 
regions becomes negative. Stars, which are preferentially on radial 
orbits, escape leaving mainly stars on tangential orbits. Finally, just 
before cluster disruption, the anisotropy in the whole system be¬ 
comes slightly positive. At that time of cluster evolution most stars 
have radial orbits. Clusters concentrated stronger (W7235), develop 
from the very beginning a small positive anisotropy in the outer and 
middle parts of the system (see Figure 6). In the course of evolution 
the amount of anisotropy in the outer parts of the system is reduced 
substantially by tidal stripping, as stars on radial orbits escape pref¬ 
erentially. As tidal stripping exposes deeper and deeper parts of 
the system, the anisotropy (for large Lagrangian radii) gradually 
decreases and eventually becomes negative. At the same time the 
anisotropy in the middle and inner parts of the system stays close 
to zero. 

It is interesting to note that for clusters which undergo core 
collapse (W5235 and W7235) the anisotropy in the inner parts of 
the system becomes slightly positive just around the collapse time. 
This behaviour can be attributed to the binary activities, which are 


Figure 4. Evolution of the anisotropy for 10%, 40% and 99% Lagrangian 
radii for model W3235. 


0.4 


0.2 

I 0 
< 

- 0.2 


-0.4 

0 5000 10000 15000 20000 25000 30000 

Time (Myr) 

Figure 5. Evolution of the anisotropy for 10%, 40% and 99% Lagrangian 
radii for model W5235. 

mainly concentrated in the cluster core. Interactions between bina¬ 
ries and field stars (relaxation and energy generation) put stars on 
more radial orbits. 

The anisotropy of main-sequence stars shows very similar 
behaviour to that discussed above. Mass segregation forces white 
dwarfs and neutron stars (most massive stars) to preferentially oc¬ 
cupy the inner parts of the system. Therefore the anisotropy for 
them is much more modest than for the main-sequence stars, and 
stays close to zero. Anisotropy evolution agrees very well with the 
results of Paper II and qualitatively with the results obtained by 
Takahashi (1997) and Takahashi & Lee (2000). 

3.4 Mass segregation. 

Three main effects may be expected to govern the evolution of the 
mass function in models studied in this paper. First, there is a pe¬ 
riod of violent mass loss due to sellar evolution of the most mas¬ 
sive stars. It takes place mainly during the first few hundred million 
years and its amplitude strongly depends on the slope of the mass 
function. Secondly, there is the process of rapid mass segregation, 
caused by two-body distant encounters (relaxation). It takes place 
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Figure 6 . Evolution of the anisotropy for 10%, 40% and 99% Lagrangian 
radii for model W7235. 

mainly during core collapse and is relatively unaffected by the pres¬ 
ence of a tidal field. Thirdly, there is the effect of the tidal field 
itself, which becomes more important after core collapse or even 
earlier (for models with a shallow mass function and with low con¬ 
centration). Because the tides preferentially remove stars from the 
outer parts of the system, which are mainly populated by low-mass 
stars (mass segregation), one can expect that the mean mass should 
increase as evolution proceeds, making allowance for stellar evolu¬ 
tion. 

The basic results are illustrated in Figure 7 for the total mean 
mass inside the 10% and 50% Lagrangian radii and the tidal ra¬ 
dius, and in Figure 8 for the mean mass of main-sequence stars 
and white dwarfs inside 10% Lagrangian radius for model W5235. 
The violent and strong mass loss due to stellar evolution is charac¬ 
terized by an initial decrease of the mean mass. The evolution of the 
most massive stars will remove a substantial amount of mass from 
the system and, consequently, greatly lower the mean mass in the 
whole system. This behaviour is visible also in the other models, 
but with one exception. The mean mass inside the 10% Lagrangian 
radius for model W7235 is increasing instead of decreasing. For 
this model, the initial concentration is high enough to force very 
quick mass segregation in the central part of the system. The most 
massive main sequence stars (not evolved yet): white dwarfs and 
neutron stars, sink into the center. For model W3235 for most of 
the time the mean mass stays nearly constant for the whole sys¬ 
tem. Just before the cluster dissolution it substantially increases. 
For all models tidal stripping connected mainly with removal of 
less massive stars forces the mean mass in the middle and outer 
parts of the system to steadily increase. The much faster increase 
of the mean mass inside the 10% Lagrangian radius (mass segre¬ 
gation) than in the outer Lagrangian radii (50% and 99%) during 
the collapse phase nearly stops at the core bounce time. The sub¬ 
sequent evolution is characterized by the slower rate of increase of 
the mean mass (close to the time of the core bounce the mean mass 
is practically constant). This is in good agreement with results ob¬ 
tained by Giersz 8 l Heggie (1997) for small A'^-body simulations 
and for Monte Carlo simulations presented in Paper IT The reason 
for the nearly constant mean mass is unclear. The slow increase of 
the mean mass can be attributed to the binary activity, which leads 
to removal of some low mass stars from the core. For all models 
the effect of tides manifests itself by a gradual increase of the mean 
mass in the halo. 



Time (Myr) 

Figure 7. Evolution of the mean mass inside 10%, 40% Lagrangian radii 
and Ttid for model W5235. 



Figure 8 . Evolution of the mean mass for 10% Lagrangian radius for model 
W5235; ms and wd means mean-sequence stars and white dwarfs, respec¬ 
tively. 

The evolution of the mean mass for main-sequence stars and 
white dwarfs for the inner 10% Lagrangian radii is shown on Fig¬ 
ure 8. When less and less massive main-sequence stars finish their 
evolution as less and less massive white dwarfs, the mean mass of 
white dwarfs decreases with time. The newly created less massive 
white dwarfs affect the mean mass of white dwarfs more strongly 
for the steeper mass function than for the shallower one. For the 
steeper mass function a smaller number of massive white dwarfs 
is created comparable to the shallower mass function. At the time 
around the core bounce the average mass of white dwarfs increases. 
This is connected with the decreasing core size and continuing 
mass segregation. Deeper in the cluster center a higher fraction 
of massive white dwarfs is present. In the post collapse phase the 
white dwarfs average mass is decreasing. Binaries start to be cre¬ 
ated mainly from the most massive stars (neutron stars and white 
dwarfs) deep in the core. They dynamically interact with other stars 
and in the course of evolution are removed from the system together 
with some massive white dwarfs. This leads to the decrease of the 
mean mass of white dwarfs. 

The mean mass of main-sequence stars initially slightly de- 
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Figure 9. Evolution of the mass ratio of main-sequence stars (ms), white 
dwarfs (wd) and neutron stars/black holes (ns) to the actual total cluster 
mass and the evolution of the ratio of the actual total mass (ms+wd+ns) to 
the initial total mass for model W5235. 



Figure 11. Evolution of the power-low index of the mass function for dif¬ 
ferent shells: R(0-30), R(30-60), R(60,rt) and for the global mass function 
for model W5235. The power-low index was computed for main-sequence 
stars in the mass range between O.IMq to O.SMq. 



Time (Myr) 


Figure 10. Evolution of the ratio of evolved stars to the actual total mass for 
different zones: R(O-O.Ol), R(0.01-0.05), R(005-0.1), R(0.1-0.3), R(0.3-0.6) 
and R(0.6-rt) for model W5235. 


creases, which is connected with stellar evolution of the most mas¬ 
sive stars. Then, there is a period of the gradual increase of the 
mean mass due to mass segregation. At late phases of cluster evo¬ 
lution the rate of increase of the mean mass speeds up. Around the 
core bounce time one can observe the break of the rate of increase 
of the main-sequence average mass. This can be again attributed to 
binary activities. 

The evolution of the total mean stellar mass in the different 
places of the system (Figure 7) and the evolution of the average 
mass inside 10% Lagrangian radius for main-sequence stars and 
white dwarfs (Figure 8) agree very well with the results obtained in 
Paper II. 

Stellar evolution will change a stellar content of the cluster. 
In the course of time main-sequence stars will evolve creating first 
black holes, neutron stars and then white dwarfs with decreasing 
masses. In Figure 9, evolution of the mass ratio of main-sequence 
stars (ms), white dwarfs (wd) and neutron stars/black holes (ns) 
to the actual total cluster mass and the evolution of the ratio of 
the actual total mass (ms + wd + ns) to the initial total mass. 


are shown for model W5235. The fraction of evolved stars present 
in the system at any time results from two processes which act in 
opposite directions: stellar evolution (formation) and evaporation 
across the tidal boundary (loss). At the beginning of cluster evolu¬ 
tion the sharp increase of the mass ratio of evolved stars and de¬ 
crease of the mass ratio of main-sequence stars and the total mass 
ratio are connected with the stellar evolution of the most massive 
stars in the system. Then, the steady increase of the mass ratio for 
wd and ns is observed. After the time of core bounce the ratio 
for ns is slowly decreasing because of binary activities (the most 
massive stars (neutron stars) are with the highest probability in¬ 
volved in binary formation and finally removed from the system by 
interactions with the field stars and other binaries). The ratios for 
ms + wd + ns and ms are decreasing, as expected. At the time 
close to the dissolution time evolved stars consists of up to 75% 
of the mass of the all stars and up to 63% of the number of all 
stars. The fraction of mass/number of wd is larger for systems with 
shorter relaxation time and with stronger influence of the tidal field 
of a parent galaxy. This is in a very good agreement with the results 
obtained by VH. 

Figure 10 shows the evolution of the mass ratio of evolved 
stars to the actual total mass for different zones in the system (see 
description of Figure 10). As expected, the mass segregation pro¬ 
cess and removal of less massive stars due to the evaporation pro¬ 
cess across the tidal boundary cause steady increase of the mass ra¬ 
tios. For the inner part of the system (up to about 5% Lagrangian ra¬ 
dius), during the core collapse time the rate of increase of the mass 
ratio is the fastest in the system. However, after the core bounce 
the rate becomes slowest The mass ratio is nearly constant. This is 
connected with the fact that the mass segregation process is nearly 
completed at the time of core bounce (see Figure 7 and its discus¬ 
sion). Further increase of the mass ratio can be mainly attributed to 
binary activities. For parts of the system outside 10% Lagrangian 
radius the behaviour of the mass ratio is opposite to that discussed 
earlier. During the collapse phase the rate of mass ratio increase is 
slower than in the post collapse phase. The further out 10% La¬ 
grangian radius the ratios are smaller and faster, respectively. This 
behaviour can be attributed to the properties of the stripping pro¬ 
cess, which is more and more effective when deeper and deeper 
parts of the system are exposed. In between these two zones (La- 
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grangian radius between 5% and 10%) the rate of increase of mass 
ratio is nearly constant during the whole evolution. The amount of 
mass in evolved stars in the core can be as high as 97% and in the 
outer parts of the system up to 55%, at the time close to the dis¬ 
solution time. The same qualitative behaviour is observed for other 
models. However, for models with steeper initial mass function and 
for more concentrated ones the mass ratios are smaller. For initially 
more concentrated models the time of core bounce is much shorter 
than for the less concentrated models. So, the process of mass seg¬ 
regation is less advanced and therefore the mass ratios are smaller. 
Also, for more concentrated cluster the evaporation process is less 
effective. The ratios are about 10% smaller than for W5325 model. 
For the model W3235, which does not enter the post-collapse evo¬ 
lution before the cluster dissolution the mass ratios behave like for 
other models in the collapse phase, but the spread between ratios 
in different zones is much smaller. It is worth to note that cluster 
at the time 15 Gyr can consist of up to 90% of evolved stars in 
the core and up to 10% of evolved stars in the outer parts of the 
system. This conclusion has very important observational conse¬ 
quences. The presence of a substantial number of practically in¬ 
visible stars has to be taken into account when the global globular 
cluster parameters are drown from the observations. 

Now lets discuss the evolution of the local (LMF) and global 
(GMF) mass functions in comparison to the initial mass function 
(IMF). According to VH the mass slope of a mass function can be 
obtained from the numerical data in the following way: 

^ _ ln{dN/dmi)i — ln{dN/dm 2)2 
ln{mi) — ln(m 2 ) 

{dN/dm)i ^2 is estimated by Afi, 2 (f)/Ami_ 2 , where A^i, 2 (f) is the 
total number of stars between mi ,2 and mi ,2 + Ami, 2 . According 
to VH nil and m 2 were chosen to O.IMq and O.5M0, respec¬ 
tively. The LMF were calculated for three zones (VH): up to 30% 
Lagrangian radius, between 30% and 60% Lagrangian radii and 
from 60% Lagrangian radius up to tidal radius. In the case when 
only mass segregation is taken into account, the LMF in the inner 
parts of the system becomes flatter and in the outer parts becomes 
steeper. The evaporation process acts in the opposite direction. It 
tends to flatten the mass function because of the preferential re¬ 
moval from the system (particularly from the outer parts) of low 
mass stars. As it can be seen in Figure 11 the theoretically predicted 
behaviour of LMF is qualitatively represented by the numerical re¬ 
sults. For all the discussed models the LMF for the outermost parts 
of the system initially increases. The increase is slightly higher for 
more concentrated systems with the same IMF (this is consistent 
with the picture of mass segregation), and even more pronounced 
for the system with the same concentration but with stepper mass 
function. For all models the GMF closely resembles the LMF in 
the middle parts of the system. The IMF during dynamical clus¬ 
ter evolution is very quickly forgotten and practically impossible to 
recover from the observational data. The actual GMF of globular 
clusters can be recovered from observational data for middle parts 
of the system (close to the half-mass radius). The results discussed 
above are in a very good agreement with the results obtained by 
VH for A^-body runs. 


4 CONCLUSIONS. 

This paper is a continuation of Paper I and II, in which it was 
shown that the Monte Carlo method is a robust scheme to study, 
in an effective way, the evolution of very large A^-body systems. 


The Monte Carlo method describes in a proper way the graininess 
of the gravitational field and the stochasticity of real A-body sys¬ 
tems. It provides, in almost as much detail as A-body simulations, 
information about the movement of any object in the system. In 
that respect the Monte Carlo scheme can be regarded as a method 
which lies between direct A-body and Fokker-Planck models and 
combines most advantages of both these methods. This is the first 
important step in the direction of simulating the evolution of a real 
globular cluster. 

It was shown that the results obtained in this paper are in qual¬ 
itative agreement with these presented by CW, VH, AH, TPZ, JNR, 
Paper II, Baumgardt & Makino (2002) and Lamers et al 2005. 
Particularly good agreement is obtained with VH’s A-body sim¬ 
ulations and what is not surprising, with results of Monte Carlo 
simulations presented in Paper IT All models survive the phase of 
rapid mass loss and then undergo core collapse and then subse¬ 
quent post-collapse expansion (except model W3235) in a man¬ 
ner similar to isolated models. The expansion phase is eventually 
reversed when tidal limitation becomes important. As in isolated 
models, mass segregation substantially slows down by the end of 
core collapse. Mass loss connected with stellar evolution dominates 
the initial phase of cluster evolution. Then the rate of stellar evo¬ 
lution substantially slows down and escape due to tidal stripping 
takes over. During this phase of evolution the rate of mass loss is 
nearly constant, and higher for shallower mass functions. Energy 
carried away by stellar evolution events dominates the energy loss 
due to tidal stripping, even though the tidal mass loss is higher. For 
the first time, because of the large number of particles in simula¬ 
tions (1, 000, 000) which results in substantial reduction of statisti¬ 
cal fluctuations, it was possible to compute in an accurate way the 
evolution of density profiles. The development of power-low den¬ 
sity profile is clearly visible and agrees very well with theoretical 
predictions. The observed power-low index is equal to about —2.2. 
The strongly concentrated model (W7235), shows a modest initial 
build up of anisotropy in the outer parts of the system. As tidal strip¬ 
ping exposes the inner parts of the system, the anisotropy gradually 
decreases and eventually becomes slightly negative. Model W5235 
from the very beginning develops in the outer parts of the system 
negative anisotropy. The cluster is not concentrated enough to pre¬ 
vent removal of stars which are preferentially on radial orbits. The 
negative anisotropy stays negative until the time of cluster disrup¬ 
tion, when it becomes slightly positive (during cluster disruption 
most stars are on radial orbits). Because of mass segregation, and 
due to evaporation across the tidal boundary, which preferentially 
removes low mass stars, the mean mass in the cluster increases with 
time. During the core collapse the rate of increase of the mean mass 
is highest in the central parts of the system (mass segregation). Af¬ 
ter the core bounce there is a substantial increase in the mean mass 
in the middle and outer parts of the system (tidal stripping), and 
more modest increase in the inner parts of the system, which is 
mainly connected with binary activity. The fraction of evolved stars 
is increasing during the cluster evolution. This fraction is larger for 
systems with shorter relaxation time and stronger influence of the 
tidal field of a parent galaxy. The mass ratio of evolved stars can be 
as high as nearly 65% at the outer parts of the system up to nearly 
98% in the core. At the time of 15 Gyr these ratios are 90% and 
10%, respectively. The presence of a substantial number of practi¬ 
cally invisible stars has very important consequences for the inter¬ 
pretation of observational data and it has to be taken into account 
when the global globular cluster parameters are drown from the 
observations. Because of stellar evolution, mass segregation and 
evaporation of stars the IMF is quickly forgotten and impossible 
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to recover from the observational data. The actual GMF closely 
resembles the LMF for the middle parts of the system (close to 
half-mass radius). These results are in an excellent agreement with 
A^-body simulations presented by VH. 

In order to perform simulations of real globular clusters the 
description of some processes, already included in the code has 
to be improved, and several additional physical processes have to 
be added to the code. Stellar and binary evolution, more accurate 
treatment for energy generation by binaries, particularly in binary¬ 
binary interaction, and proper treatment of the escape process in 
the presence of a tidal field are still waiting for improvement. One 
of the population synthesis codes developed by Hurley (Hurley et 
al 2000, 2002), Portegies Zwart (Portegies Zwart & Verbunt 1966) 
and Belczynski (Belczyhski, Kalogera & Bulik 2002) can be used 
to follow more accurately single star and binary evolution. The im¬ 
plementation of techniques used in Aarseth’s NBODY codes for 
direct integrations of a few body subsystems can cure the problem 
of energy generation in 3-body and 4-body interactions (Giersz & 
Spurzem 2003). The treatment of escapers proposed by Spurzem et 
al (2005) can be implemented to solve the problem of escapers in 
the tidal field (problem which is inherently connected with gaseous, 
Fokker-Planck and Monte Carlo codes). The tidal shock heating 
of the cluster due to passages through the Galactic disk, interac¬ 
tion with the bulge, shock-induced relaxation, primordial binaries, 
physical collisions between single stars and binaries are some of 
the processes, which are waiting for implementation into the code. 
The inclusion of all these processes does not pose a fundamental 
theoretical challenge, but is rather complicated from the technical 
point of view. The international collaboration called MODEST was 
setup a few years ago to solve problems with merging all available 
codes (hydrodynamical, stelar evolution and dynamical) into a code 
which can deal in detail with the evolution of real globular clus¬ 
ters (see http://www.iiianybody.org/inodesti. Inclusion into the 
Monte Carlo code of as much as possible physical processes will 
allow to perform detailed comparison between simulations and ob¬ 
served properties of globular clusters, and will also help to under¬ 
stand the conditions of globular cluster formation and explain how 
peculiar objects observed in clusters can be formed. These types 
of simulations will also help us to introduce, in a proper way, into 
future A-body simulations all the necessary processes required to 
simulate the evolution of real globular clusters on a star-by-star 
basis from their birth to their death. 
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